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Abstract 

We consider cosmologies in which a dark-energy scalar field interacts with cold dark matter. The 
growth of perturbations is followed beyond the linear level by means of the time-renormalization- 
group method, which is extended to describe a multi-component matter sector. Even in the absence 
of the extra interaction, a scale-dependent bias is generated as a consequence of the different initial 
conditions for baryons and dark matter after decoupling. The effect is enhanced significantly by 
the extra coupling and can be at the 2-3 percent level in the range of scales of baryonic acoustic 
oscillations. We compare our results with N-body simulations, finding very good agreement. 
PACS numbers: 95.35.+d, 95.36.+X, 98.80.Cq 
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1 Introduction 



The power spectrum of matter perturbations reflects the evolution of the Universe since the time 
of matter-radiation equality. For given initial conditions, determined by the primordial spectrum 
(usually assumed to be scale invariant), the growth of perturbations depends on the cosmological 
scenario. The calculation of the present matter power spectrum can constrain this scenario through 
the comparison of the deduced spectrum with the observed large-scale structure. A major technical 
difficulty in the realization of such a program is the failure of linear perturbation theory to describe 
present-day fluctuations with characteristic length scales below roughly fOO Mpc. At length scales 
below about 10 Mpc, the evolution is highly non-linear, so that only numerical N-body simulations 
can capture the dynamics of the formation of galaxies and clusters of galaxies. Fluctuations with 
length scales of around 100 Mpc fall within the mildly non-linear regime, for which analytical methods 
have been developed. These scales (corresponding to wavenumbers in the 0.03 — 0.25/i/Mpc range) 
are of particular interest, because they correspond to the sound horizon at decoupling, which can be 
determined by reconstructing the oscillatory behavior of the matter power spectrum due to baryonic 
acoustic oscillations (BAO). 

The various analytical methods [TJ-[2] essentially amount to resummations of subsets of pertur- 
bative diagrams of arbitrarily high order, in a way analogous to the renormalization group (RG). In 
this work we shall follow the approach of [5], named time-RG or TRG. In the context of the RG 
the various observables depend on a characteristic energy scale, and evolve as this scale is varied. 
The TRG uses time as the flow parameter that describes the evolution of physical quantities, such 
as the spectrum of perturbations. The method is characterized by conceptual simplicity. It has been 
applied to ACDM and quintessence cosmologies \E\, as well as models with massive neutrinos [BJ- 

The fundamental equations in the TRG approach are the "equations of motion", i.e. the conti- 
nuity, Euler and Poisson equations. From these, equations can be derived for the time evolution of 
correlation functions for the density and velocity fields. The various spectra are obtained through 
appropriate Fourier transforms of the correlation functions. The method results in a coupled infinite 
system of integro-differential equations for the time evolution of the spectrum, bispectrum etc. The 
crucial approximation, that makes a solution possible, is to neglect the effect of the higher-order 
correlation functions in the evolution equations of the lower-order ones. The calculations performed 
so far take into account the spectrum and bispectrum and set the higher-level spectra to zero. 

The procedure of truncating the system of equations is commonly employed in the applications of 
the Wilsonian RG to field theory or statistical physics. (For a review, see [11] .1 The accuracy of the 
calculation can be determined either by enlarging the truncated system (by including the trispectrum, 
for example) and examining the stability of the results, or by comparing with alternative methods. 
The second approach is often followed, because enlarging the truncation can increase the complexity 
of the calculation considerably. In the case of the TRG, the agreement with results from N-body 
simulations for ACDM has been confirmed [5]. Also, a comparative analysis of several analytical 
methods, using N-body simulations as a reference, has been carried out in ref. [10] . The study 
demonstrates that the TRG remains accurate at the 1-2% level over the whole BAO range at all 
redshifts. 

In fig. 5 of ref. [10] the deviation of the TRG prediction from a reference spectrum derived 
through simulations for ACDM is depicted. At a redshift z = 1 the deviation is at the 1% level, or 
smaller, over the whole depicted range < k < 0.15/i/Mpc. For z — the deviation may exceed 
1% for some values of fc, but stays below 2% over the whole depicted range. Based on these findings 
and the analysis of [5] we estimate an accuracy of 1% for z = 1 and 2% for z = over the range 
< k < 0.2 /i/Mpc. The accuracy of 1-loop standard perturbation theory (SPT) can be inferred from 
fig. 1 of ref. [10] , At z — 1 the accuracy is at the 1% level for < k < 0.1 /i/Mpc, while at z = it 
is at the 2% level for < k < 0.05 fo/Mpc. We depict these ranges in figs. [3HpJof the current paper. 
The additional approximations that we make in this paper for the study of models with non-zero 
coupling between dark matter and dark energy induce uncertainties at the sub-percent level. We 
discuss this issue in detail in subsection 3.1. For this reason the level of accuracy of our results for 
the power spectra in the coupled case is expected to be similar to that for ACDM. Its magnitude is 
set by the truncations in the evolution equations (the omission of the effect of the trispectrum and 
higher spectra), which are similar in all models. 

The purpose of the present work is threefold: 
1) We extend the formalism to more complicated models. We introduce two modificiations to previous 
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studies: a) Within the matter sector we allow for an arbitrary number of species with independent 
spectra. These include baryonic matter (BM) and cold (pressureless) dark matter (CDM). One may 
also consider contributions from massive neutrinos etc. b) We allow for an interaction between CDM 
and dark energy (DE). We consider a class of quintessence models, in which there is direct coupling 
between CDM and the quintessence field. The form of the interaction is a generalization of the 
universal coupling to all species present in scalar-tensor theories in the Einstein frame. It is modelled 
through the dependence of the mass of the CDM particles on the quintessence field |12j-|20], 

2) We test the accuracy of the method in this enlarged framework by comparing with available N- 
body simulations. We perform our numerical analysis for a model for which results from simulations 
are given in ref. [3D]. In the context of coupled quintessence, the cosmological evolution can be very 
diverse 14, 201 1211 [22] , It is very time-consuming to study exhaustively every model through N-body 
simulations. Our approach provides an alternative method, which can be much faster, while retaining 
the necessary accuracy. It is also important to note that, while the N-body simulations are highly 
accurate at the length scales of galaxies and clusters of galaxies, they are less accurate in the BAO 
range because of the required large volumes. On the other hand, analytical methods, like ours, are 
more accurate in the quasi-linear regime of large length scales. The two approaches can be viewed as 
complementary. 

3) We provide predictions for observables for which non-linear corrections can be important. As 
such, we study the bias between dark and baryonic matter in the BAO range for models of coupled 
quintessence. 

The couplings between the matter sector and DE are constrained by observations. For the BM- 
DE coupling, the bound from the Cassini spacecraft 23 limits its order of magnitude to be below 
10 -3 . As such small couplings produce negligible effects on the power spectrum, we assume that 
the BM-DE coupling is exactly zero. The coupling between CDM and DE is constrained by various 
considerations, such as the modification of the spectrum of the cosmic microwave background (CMB) 
or that of the matter distribution. A common feature of the class of models we are considering is 
that the presence of an additional long-range force between CDM particles, induced by the DE field, 
modifies their clustering properties. The various observable consequences have been discussed in the 
literature [24]-[40]. The strength of the CDM-DE interaction is constrained through the comparison 
with the observed CMB and matter spectra. It has been shown that, for particular models, the CDM- 
DE coupling must be considerably smaller than the gravitational one [381139] . Such constraints cannot 
be considered generic, because the evolution of the cosmological background and the perturbations 
around it varies considerably from model to model. We work within the model of [40| . because our 
main objective is to compare with the results of N-body simulations presented there. The couplings 
that we consider are roughly consistent with the bounds deduced in [381 139] for a model similar to 
ours. 

In this work we follow a novel approach for the derivation of the fundamental equations. We 
derive the continuity, Euler and Poisson equations on an expanding background, starting from the 
conservation of the energy-momentum tensor. In order to cast these equations in a form that general- 
izes the standard expressions on a static background, we need to impose a certain hierarchy between 
the density perturbations, the velocity field and the potentials. Our derivation makes it straightfor- 
ward to generalize these equations in future studies in order to take into account non-zero pressure 
and higher-order terms. Next, we derive the system of differential equations for the spectrum and 
bispectrum, within a truncation that neglects higher-level spectra. We integrate these equations nu- 
merically in order to produce the non-linear spectra at low redshift and compare with the results of 
N-body simulations. We study in detail the difference, usually characterized as bias, between the 
spectra of dark and baryonic matter. We find that the CDM-DE coupling enhances significantly the 
bias of the decoupled case. 

In the following section we derive the evolution equations for the spectra of dark and baryonic 
matter. The details of the derivation for the case of one massive component are given in appendix 
A. The generalization for an arbitrary number of massive components is presented in appendix B. 
The results of the numerical integration of the evolution equations are presented in section 3. We 
compare them with the results of N-body simulation. We also discuss in detail the form of the bias 
in the BAO range. 
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2 Dark matter coupled to dark energy 



2.1 Non-linear evolution equations for the perturbations 

We assume that the energy density of the Universe receives significant contributions from three 
components: a) standard baryonic matter (BM); b) a species of weakly interacting, massive particles, 
which we identify with cold dark matter (CDM); and c) a slowly varying, classical scalar field <f>, 
whose contribution to the energy density is characterized as dark energy (DE). We also consider the 
possibility that there is a direct coupling between the CDM particles and the scalar field. Its equation 
of motion takes the form 



1 8 ( , — M „ 84) \ dU 



g dx^ V dx v 



--TT+m {Todm)*. (2.1) 



We normalize all dimensionful quantities, such as the scalar field, with respect to the reduced Planck 
mass M = (87rG) _1//2 . The full M-dependence is displayed explicitly in appendix A. Our normaliza- 
tion here is equivalent to setting M = 1. Equation (|2.1[) can be obtained if we assume that the mass 
m of the particles has a dependence on <f> [JT]. Then we have P((f>) = — dhxm(cj>)/d(j). In order to be 
consistent with the stringent observational constraints for the baryonic sector, we assume that the 
interaction with the DE scalar field is confined to the CDM sector. The BM has no direct coupling 
to (j>. 

For the metric, we consider an ansatz of the form 

ds 2 = a 2 (r) [(1 + 2$(r,f))dr 2 - (1 - 2*(t, x)) dx dx] . (2.2) 

We assume that the Newtonian potential $ is weak, $ <C 1, and that the field 4> can De decomposed 
as 

4>(T,x) = $(T) + 5<t>(r,x), (2.3) 

with d<f>/<t> <C 1. In general, <f> = 0(1) in units of M. The magnitude of the fluctuations of <j> is 
expected to be similar to that of the gravitational field The reason is that the source for both is 
the dark matter density, to which they couple with comparable strength (as will be apparent in the 
following). Finally, the density can be decomposed as 

p(r,x) = p{t) + 5 P (t,x). (2.4) 

We allow for significant density fluctuations, even though our analysis is not applicable when they 
are much larger than the background density. Our aim is to take into account the effect of the local 
velocity field Sv, when this becomes significant because of large field gradients. For subhorizon per- 
turbations with momenta k ^> H = a/a, the linear analysis predicts \8v\ ~ (k/l-L)<& ~ (H/k)(8p/ p). 
A consistent expansion scheme can be obtained if we assume that <3> \Sv\ <S 1. Including the den- 
sity perturbations, our assumptions can be summarized in the hierarchy of scales: $, 5<f>/<j> <C \5v\ 
Sp/p < 1. At the linear level, we have Sv 2 ~ &(6p/p). We assume that such a relation holds at the 
non-linear level as well, within the range of applicability of our scheme. The velocity field is driven 
by the spatial derivatives of the potentials $, 84>/<j>. As we are dealing with subhorizon perturbations, 
it is consistent to make the additional assumption that the spatial derivatives of 5>, 8(j> dominate over 
their time derivatives. The predictions of the linear analysis allow us to make a more quantitative 
statement. We assume that a spatial derivative acting on 8<j> or Sv increases the position of that 
quantity in the hierarchy by one level. In this sense V"I> is comparable to 8v, while V 2< 3? is comparable 
to p. 

With the above assumptions, one can derive the equations that describe the evolution of the 
Universe. For one non-relativistic species, the derivation is presented in appendix A. It is generalized 
to two species in appendix B. The evolution of the homogeneous background is described by 

H 2 =^[a 2 P* + \^ + a2[/ (^)] = ^a 2 p to t (2.5a) 

i=l,2 

pi + 3H P i = -fr $ Pi (2.5b) 
$ + 2H$ = -a 2 E > ( 2 - 5c ) 
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where we have defined p tot = J2 t Pi + tfl&a 2 ) + U(<f>). For the CDM we have a non-zero constant 
Pi = P, while for BM, because of the strong constraints from solar system tests of General Relativity 
[23] . we set /3 2 = 0. 

We describe the perturbations in terms of the scalar-field perturbation 8(f), the Newtonian potential 
<&, the density perturbations 5 pi and the velocity fields Vi. We have two Poisson equations 



X7 2 5</> = - a 2 frSpi = PiW'ttiti 

i i 



(2.6a) 
(2.6b) 



with Qj(r) = pia 2 /(3H 2 ), and the continuity and Euler equations 



5p\ + ZHSpi + V[(p; + 8pi)5vi] = — Pa, 



(2.7a) 
(2.7b) 



2.2 The CDM-BM quadruplet and the power spectra 

The evolution equations are expressed in their most useful form in terms of the density contrasts 
5i = 5 pi I pi < 1 and 0i(k,r) = V ■ 5vi(k,r). For the Fourier transformed quantities, we obtain from 
eq. (|2.7a[) 



6i(k, t) + BiQt, t) + J d 3 k! d 3 k 2 S D {k - ki - k 2 ) a(ki, ka) 4 (k x , r) <5i(k 2) r) = 0, 



where 



a(ki, k 2 ) = 



ki ■ (ki + k 2 



Eqs. (f277b|) . (|2T6a|) . (|2T6b)l give 



0j(k,r)+CH-ft(/>)0i(k,T) + 



7" 



+ / d 3 k 1 d 3 k 2 5 I5 (k-k 1 -k2)/3(k 1 ,k2)6l l (k 1 ,r)^(k 2 ,r) = 0, 



where 



/3(ki,k 2 ) = 



(ki + kafkx ■ k 2 



(2.8) 
(2.9) 

(2.10) 
(2.11) 



In appendix B we discuss the above equations for an arbitrary number of non-relativistic species. 
Here we concentrate on the case of interest, i.e. CDM coupled to </> and BM only gravitationally 
coupled. We define the quadruplet 



/ Vi(k,»?) \ 
V2(k,r?) 
<£3(k,?7) 

V <^4(k,77) / 



6cDM(M,r]) 
U 

5 BA f(k, rf) 
6> BA /(k, rj) 



\ 



(2.12) 



H 



J 



where r\ = In a(r). This allows us to bring eqs. (|2.8|) . (|2.10|) in the form [HUE] 

d v (fia0<.,rj) + £lab<fib(k., -q) = e r ''y abc (k, -ki, -k 2 )<£i,(ki, »7)<£ c (k 2 , rf). 



(2.13) 



The indices a, b, c take values 1, . . . , 4. The values 1,2 characterize CDM density and velocity per- 
turbations, while 3,4 refer to BM quantities. Repeated momenta are integrated over, while repeated 
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indices are summed over. The functions 7, that determine effective vertices, are analogous to those 
employed in [4J [5] . The non-zero components are 



7121 (k, ki, k 2 



5(ki,ka] 



5 D (k + ki +k 2 ) = 7n 2 (k,k 2 ,ki) 



The fi-matrix is 



7 222 (k,ki,k 2 ) = /3(ki,k 2 ) 5 D (k + ki + k 2 ) 

7343(k, k 3 ,k 4 ) = " ' Sr> (k + k 3 + k 4 ) = 7 334 (k, k 4 , k 3 ) 



744 4 (k, k 3 , k 4 ) = /3(k 3 ,k 4 ) 5 D (k + k 3 + k 4 ) 
/ 1 - 1 



(2.14) 



"(»?) = 







ft 



V 







-1 



(2.15) 



(2.16) 



where a prime denotes a derivative with respect to 77. 

The next step is to derive evolution equations for the power spectra. The spectra, bispectra and 
trispectra are defined as 

(Va(K v)<Pb(ci, 77)) =5 D (k + q)P ai ,(k, 77) 
(tpaQt, v)<Pb(<l, il)<p c (p, 7j)> =<5z>(k + q + p)Babc(k, q, p, 77) 
(ip a (K rj)<Pb(*l, r?)yc(p, r?)pd(r, tj)> =fo(k + q)<5 D (p + r)P ab (k, r))P c d(p, v) 

+ 5 D (k + p)5 D (q + r)P ac (k, J7)P M (q, 77) 
+ (k + r)5n (q + p)P a<i (k, 77)P 6c (q, 77) 
+ S D (k + p + q + r)Q abcd (k, p, q, r, 77) . 

In appendix B we summarize the derivation of the evolution equations for an arbitrary number of 
species. The essential approximation that we have to make in order to obtain a closed system is to 
neglect the effect of the trispectrum on the evolution of the bispectrum. In this way we obtain 

d v P a b(k,rj) = -Q ac P cb (k, 77) - D. bc P ac (k, 77) 

+e v J d 3 q[y acd (k, q, q - k)B bcd (k, q, q - k) 

+Jbcd{K q, q - k)B acd (k, q, q - k)] , (2.17) 
-~tt ad B dbc (k, q, q k) - Q bd B adc (k, -q, q- k) - Q cd B abd (k, -q, q - k) 

+2e" f d 3 g[ 7a(i(! (k,-q,q-k)P d i,(q,77)P !C (k-q,77) 



(2.18) 



d v B abc (k,-q,q- k) 



+Jbde(— q, q- k,k)P dc (k- q, 77)P ea (k, 77) 
+7cde(q- k,k, -q)P da (k,77)P e6 (q,77)l. 
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The formal solution of the above system is given by 



£f ac (k , 77, 0) g bd (k , 77, 0)P cd (k , r? = 0) 

+ / I d3q9ae ( k ' v, ' ri ' S>9bf ( k, ' ri, ' ri ' S> 

x [7<:cd(k, q, q - k) B fcd (k, q, q k; 77') 

+ 7/cd(k, q, q - k) B ecd (k, q, q k; 7/)] , (2.19) 

g ad (k,r),0)g be (-q,r),0)g cf (<i-k,ri,0)B <ie f(k, q, q k; 77 = 0) 
/•1 , 

+2 / dri'e** g ad (k , 77, n')9be{-q. , V, v')9cf(q - k , r), n ) 
Jo 

x [ld 9 h(k, q, q- k)P e g(q,ri')Pf h (q- k,r)') 
+le 3 h{-q., q k, k)P/ g (q - k , 7/)P dh (k , 77') 
+ 7/gh (q - k, k, -q)Pd 9 (k , V ')P,h((l , V')] , (2.20) 

where g a b(k , 77, 77') is the linear propagator, which gives the evolution of the field at the linear level: 
<pi(k,n) =g ab (k,i 1 ,r ] ')^(k,i 1 '). 

The solutions can be expanded in powers of the interaction vertex j a bc, in order to establish the 
connection with perturbation theory [5]. The lowest order, corresponding to linear theory, is obtained 
by setting j a t c — 0. The linear spectrum P^ b and bispectrum B^ bc are given by the first line of each of 
the above equations. The 0(7) correction for the bispectrum is obtained by inserting P^ b in place of 
P a b in the r.h.s. of eq. (|2.20p . Inserting the bispectrum at this order in eq. (|2.19|l generates the 0(7) 
and 0(-y 2 ) contributions to the power spectrum. At this order, the result for the power spectrum 
reproduces exactly the one- loop expression in standard perturbation theory (SPT) [42]. Iterating 
the procedure generates the higher-order corrections. However, differences with perturbation theory 
arise at higher orders, because of the approximation Q a bcd = that we have made in deriving the 
evolution equations for the power spectrum and bispectrum. 

3 Numerical analysis 
3.1 Approximations 

The presence of two massive species (BM and CDM) complicates the structure of the equations 
compared to the case where they are treated as a single fluid, discussed in |5]. The full system of 
eqs. (|2.17|l . (I2.18|) contains 74 equations, namely, 10 for the power spectra and 64 for the bispectra, 
compared to the 11 equations of the single-matter case. An accurate calculation also requires the 
discretization of the fc-space with at least 500 points. Taking into account that the bispectra depend 
on three external momenta, it is apparent that the necessary computing power is significant. 

The system can be reduced if additional approximations are made, based on the following obser- 
vations: 

• The dynamical vertices of eq. (|2.14fl do not mix the CDM components with the BM ones. The 
coupling between the two type of components is entirely due to the linear part of the equations, 
and especially to the O23 and Q41 entries of eq. <|2. 15|) . through which the 1,2 and 3,4 indices 
are mixed. These originate in the Poisson equation, in which the fluctuations of all the matter 
species contribute universally to the gravitational potential. 

• The ratio of BM and CDM density perturbations is usually characterized as bias: b = 8bm /Scdm- 
If at early times b is independent of k, the subsequent linear evolution preserves this indepen- 
dence, so that b is only a function of 77. At the linear level, the density-velocity and velocity- 
velocity spectra are proportional to the density-density ones, with fc-independent proportionality 
factors. These factors are appropriate powers of the linear growth functions, so that the growing 
modes are selected [5]- 

• The effect of the CDM-DE coupling /3 on the bias can be estimated analytically in the linear 
approximation. For /3 = the evolution equations for CDM and BM are the same, so that, for 
identical initial conditions, we have Pn = P33 = P13. For /3 7^ 0, the growing mode solution 



Pab (k, 77) = 



B abc (k, -q, q - k; 77) = 
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of the linearized equations has the form ip a (k,r]) = <p(k, rf) [1, f(il),b(ri),b(ri)f(ri)], where /(??), 
b{rf) are solutions of the system 

^Q C Dm{2P 2 + 1) + |n SM 6- \\^CDM - ~P,BM +P4>'f = (3.1) 

f + J[f + f + f 2 -\\^CDM-\n BM = 0, (3.2) 

with the prime denoting derivatives with respect to 77 = In a(r). For f3 = 0, we have b = 60 = 1 
and / = foijf), with fo(v) the growth function of the corresponding decoupled model. The 
corrections for f3 7^ are 0(/3 2 ). This is obvious for those arising from the first term in eq. 
(|3.1[) . The last term in the same equation has a similar effect because the evolution of the field 
is given by eq. (|2.5c|l . In all models in which the CDM-DE coupling affects the cosmological 
evolution, the two terms in the r.h.s. of eq. (|2.5c[) are comparable, so that <j>' = O(P). 

• The non-linear corrections induce a fc-dependence at low redshifts and large k. Again this is 
an 0{/3 2 ) effect. This can be verified in the context of the loop expansion, as the vertices are 
/3-independent, while the propagators for CDM and BM differ by terms of 0(/3 2 ). 

• Finally, there is a bias induced by the initial conditions for the CDM and BM spectra, which 
are not identical at the end of the decoupling era. The magnitude of this effect can be deduced 
from our subsequent analysis, and is apparent in fig. HJ in which the bias in the BAO region 
is depicted at a redshift z = 1.12. The bias deviates from 1 at a level smaller than 2% for the 
decoupled scenario (fi = 0). We have checked that at z = the effect is below 1%. 

• The estimates on the /3-dependence of the bias are also confirmed by fig. [6] Within the BAO 
range, the bias factor b receives corrections of 0(/3 2 ) relative to the /3 = case. 

These observations offer the possibility to reduce the number of evolution equations by computing 
only the pure CDM or BM spectra and approximating the mixed ones. In the resulting equations for 
the pure spectra we make the following approximations: 

• In the r.h.s. of eqs. (|2.17|l we approximate the mixed power spectra as the geometrical averages 
of the pure ones, i.e. 

Pis(M) ^ s/Pii(k,r))P 3 3(k,r]), Pu(k,ri) ~ s/Pii(k,r))P M (k,ri), 

P 23 (fc, 77) ^ \/ p 22(k,v)P33{k,v), p 24(k, r?) ~ y/P 2 2(k,r ] )P i 4k,v) ■ (3-3) 

If the bias were fc-independent, we would have Pi 3 = iPn, P33 = b 2 Pn, and the first relation 
would be exact. The same is true for the other relations, because of the fc-independence of the 
linear growth functions. The accuracy of the resulting mixed power spectra can be estimated 
through the fc-dependence of the bias. As we have discussed above, this is an 0{/3 2 ) effect, 
which is smaller than 1% for the values of /3 that we use (/3 < 0.1). 

• In the r.h.s. of eqs. (|2.18[) the mixed components of the bispectra are approximated by the 
corresponding pure ones, e.g. 

Bus ^ Pin , (3.4) 

and so on. The direction of the approximation (whether a CDM index is turned into a BM one, 
or vice versa) is decided by a majority criterium: if a bispectrum has two CDM components and 
a BM one, it is approximated by a purely CDM bispectrum, and vice versa. The accuracy of 
this approximation is determined by the magnitude of the bias. The deviation of the bias from 
1 receives a correction around 1-2% at low redshifts because of the different initial conditions 
for CDM and BM, and a correction of 0(/3 2 ) because of the CDM-BM coupling /3 ^ 0. 
We emphasize that this approximation affects the calculation of the power spectra only indi- 
rectly. The approximated mixed bispectra appear in the r.h.s. of eq. Q2.18P along with several 
pure ones. This induces an error in the pure bispectra, obtained through the integration of eq. 
(|2.18|l . which is significantly smaller than 1-2%. The pure bispectra then affect the calculation 
of the power spectra by appearing in the r.h.s. of eq. (|2. 1 7f) . We expect the residual effect of 
approximations such as (|3.4[) on the accuracy of the computed power spectra to be below 1%. 
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Density-Density Power Spectra at z=0 ((3=0.05) 
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Figure 1: Comparison of results from N-body simulations and our calculation (/3 = 0.05). We display the 
ratio of the non- linear and linear spectra for z = 0. 

Through the above approximations the system is reduced to a set of 22 coupled equations, which 
can be solved in a way analogous to that described in Appendix B of [5]. As we have explained in 
detail, our approximations are expected to be valid at the sub-percent level for the power spectra 
within the BAO range. This can be verified a posteriori, by comparing with N-body simulations, as 
we will do in the following section. Of course, the most unambiguous test would be the comparison 
with an exact solution of the full 74 equations (|2.17[) . (|2.18|l . Unfortunately, it is very difficult to 
keep the numerical accuracy of the solution of such a huge system of integro-differential equations at 
the sub-percent level, in order to carry out this test. 

3.2 Results 

We have applied the formalism of the TRG to a quintessence model with non-zero coupling between 
CDM and the quintessence field 0. We chose a particular model for which the matter spectrum has 
been calculated in [40] through numerical N-body simulations. The field has a potential V(4>) ~ 
with a — 0.143. The present-day energy content of the Universe has S!de = 0.743, Qcdm = 0.213, 
Q.bm = 0.044. The Universe is assumed to have vanishing spatial curvature (fife = 0), current 
expansion rate Ho = 71.9 km s _1 Mpc -1 . The mass variance is taken as = 0.769, as calculated from 
the linear spectrum. 

The initial conditions for the integration of the evolution equations for the spectra have been set 
at a redshift z — 40. At such early times the evolution is linear to a very good approximation. We 
employed the implementation of the background and linear-perturbation equations in the Boltzmann 
code CMBEASY [33], generalized for the interacting case |40] , We assumed that the primordial 
spectrum is scale invariant with spectral index n = 0.963. We chose the initial value of the scalar 
field close to its tracker value in the uncoupled case, and adjusted the value of the dimensionful 
constant in its potential so as to obtain the present-day energy content listed above. 

In the ACDM case the momentum dependence of the power spectrum ensures that the momentum 
integrations in the evolution equations (|2.17p . (|2.18p are both infrared (IR) and ultraviolet (UV) 
finite at any order in perturbation theory. However, in order to perform numerical computations 
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Density-Density Power Spectra at z=0 (P=0.1) 
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Figure 2: Comparison of results from N-body simulations and our calculation (/3 = 0.1). We display the 
ratio of the non- linear and linear spectra for z = 0. 

these integrations must be cut off both in the IR and the UV. An appropriate IR limit eliminates 
contributions from very large length scales, of the order of the horizon distance. In our numerical 
study of the case of coupled quintessence we employ an IR cutoff km ~ 10 _2 /i/Mpc. We have 
checked that an IR cutoff km ~ 10 _3 /i/Mpc does not alter our results for the spectra. It tends, 
however, to increase the noise for numerical integrations with larger time steps than the ones we 
employ. In principle, the UV cutoff must also be chosen with care. Contributions with very large 
momenta correspond to length scales for which our method is not accurate. At scales below a few 
Mpc the process of virialization is crucial for the formation of galaxies and clusters of galaxies. 
The non-linear corrections that we take into account through the TRG are not sufficient for the 
quantitative description of the physics at these scales. For this reason we implement an UV cutoff 
kuv — 2.3/i/Mpc. We have checked that the variation of kuv by a factor of two induces a variation 
of the spectra in the BAO range (0.03/i/Mpc & k & 0.25/i/Mpc) at the sub-percent level. As a 
result, the evolution at small length scales does not influence appreciably the evolution at the scales 
relevant for the BAO, which are the main focus of our calculation. 

We can check the reliability of the TRG approach by comparing our results with those from 
N-body simulations. In figs. [1] and [2] we display the ratio of non-linear to linear power spectra at 
2 = 0, for P — 0.05 and 0.1 respectively. The lines correspond to our results for the CDM and BM 
density-density power spectra. The points indicate the results of the N-body simulations presented in 
[40] for the same quantities. No error bars are given for these results that would permit an accurate 
assessment of the level of agreement with our findings. On the other hand, we observe consistency of 
the two methods at scales below 0.5 ft/Mpc: The growth of the non-linear power spectrum relative 
to the linear one at small length scales is captured well by the TRG. It is also obvious from the same 
figures that the two approaches give results that start to deviate for k 0.5 ft/Mpc. This is expected, 
as the TRG cannot capture the processes of formation of bound structures that dominate at large 
momentum scales. Our main interest lies in the region 0.03/i/Mpc ^ k JS 0.25/i/Mpc, in which the 
BAO are visible. We also point out that the N-body simulations are rather noisy at such length scales 
because of finite volume effects, as they are optimized for scales smaller than the BAO range. For this 
reason, the two methods, N-body simulations and TRG, can be viewed as complementary: N-body 
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BAO Dark Matter Density-Density Power Spectra at z = for different P 
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Figure 3: Dark matter density-density spectra for various /3 at z = 0, normalized with respect to the 
smooth function of |44j . The spectra have been multiplied by an additional /3-dependent factor in order 
to be equal to 1 for k — > 0. We indicate the maximum k for which a certain level of accuracy is achieved 
by TRG and SPT. 



simulations give a reliable description of the process of virialization at relatively small length scales, 
while the TRG is the appropriate tool for the study of non-linear effects at the BAO range. The 
agreement of the two methods in the intermediate range of overlap is a confirmation their consistency. 

In figs. [3] f3] [S] we display the linear and non-linear power spectra for models with /3 — 0, 0.05 and 
0.1. We focus on the BAO range. In order to provide a clear depiction of the oscillatory behavior of 
the spectra we divide them by the smooth function of ref. [44] for £Ibm = 0. 

In fig. [3] we display the density-density power spectra for CDM for various values of f3. It 
is apparent that the inclusion of the non-linear effects causes the spectrum to increase at large 
momenta. For k 5s 0.06/i/Mpc (a region that includes the first peak), the linear and non-linear 
spectra are indistinguishable. This indicates that the non-linear corrections become negligible at 
large length scales. On the other hand, the non-linear spectrum becomes significantly larger than 
the linear one at smaller length scales. The non-linear effects are quantitatively important in the 
BAO range. In the vicinity of the second peak they provide corrections at the level of a few percent, 
while in the vicinity of the third peak the corrections are about 10%. These features of the non-linear 
corrections on the matter power spectrum in the (3 = case have been discussed in full detail in the 
past (for instance, in [5]). 

The presence of an extra CDM-DE coupling causes an additional enhancement of the spectrum. 
For example, for k — 0.3/i/Mpc the differences of the non-linear spectra for /3 = 0, 0.05 and 0.1 
have roughly doubled compared to the corresponding differences of the linear spectra. A qualitative 
change in the spectrum is that the third peak ceases to exist for ft ~ 0.1. It is likely that the 
disappearance of higher-order peaks is common in models of coupled quintessence. However, it is not 
possible to investigate this feature in a model-independent way. The large variability of the evolution 
of the cosmological background and the perturbations in models of coupled quintessence makes it 
very difficult to identify generic features. 

In fig. f?]we display the density-density power spectra of BM. They are very similar to those of 
CDM. A close inspection reveals that the enhancement of the spectrum with increasing /3 is smaller 
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gure 4: Baryonic matter density-density spectra for various /3 at z = 0, normalized similarly to those 

fig.H 

than for the CDM case. This is more apparent for k — 0.3/i/Mpc. The reason can be traced to 
the additional attractive force between CDM particles, mediated by the quintessence field. As the 
additional force is not felt by the baryons, the enhancement affects CDM and baryons differently, i.e. 
a bias is produced. This is apparent already at the linear level, but it becomes an even stronger effect 
at the non-linear level. The baryons are not subject to this force and tend to collapse more slowly 
[25] . The bias is scale dependent, as we discuss in detail below. 

In fig. [5] we depict the velocity- velocity power spectra for CDM. Similarly to the density-density 
spectra, the higher-order peaks are washed away by the combined effect of the non-linear corrections 
and the CDM-DE interaction. The non-linear density-density (dd) , density- velocity (dv) and velocity- 
velocity (vv) power spectra display similar behavior. The BM spectra follow closely the variation 
of the corresponding DM ones. There is also a hierarchy in the magnitude of the spectra for large 
k. The reason is that the velocity field is much smaller than the density perturbation at subhorizon 
scales. As a result, the non-linear effects that enhance the spectra at large k are more pronounced 
for the density-density spectra and smallest for the velocity-velocity ones. 

The bias in the BAO region is depicted in fig. [6] for a redshifts z = 1.12. It is defined as 
b(k) = (P B M(k)/P C DM(k)) 1/2 . The small kinks in the curves give an indication of the numerical 
errors in our calculation. We notice that, even at vanishing coupling (/3 = 0), there is a bias as a 
result of the different initial conditions between BM and CDM after decoupling 25 . However, in 
the /3 = case the bias parameter tends to unity at lower redshifts, as is well known. When /3 ^ 
the initial unbalance between baryons and CDM is never washed out, as a consequence of the extra 
scalar force acting on the latter and not on the former. In linear perturbation theory, this causes a 
scale-independent bias in recent epochs [25]. When non-linear corrections are taken into account, the 
situation changes dramatically. The CDM and BM components, having different initial conditions, 
evolve differently even in the /3 — case, with the BM fluid being always more "linear" than the 
CDM one. Moreover, since the non-linear growth factors are scale-dependent, non-linear corrections 
cause a scale-dependent bias even at /3 = [26] . It is clear from fig. [5] that both effects, i.e. the bias 
and its scale dependence, are enhanced by a non-vanishing coupling. This feature could provide a 
distinctive signature for this kind of models. These results are consistent with the conclusions of refs. 
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BAO Dark Matter Velocity-Velocity Power Spectra at z = for different P 
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Figure 5: Dark matter velocity-velocity spectra for various /3 at z — 0, normalized similarly to those in 

fig. m 

[401 145] . in which the emphasis is put on the halo region of the spectrum. 

4 Conclusions 

In this paper we have extended the TRG formalism introduced in ref. [5] in two respects. Firstly, we 
described the matter sector of the theory keeping track of the CDM and BM components separately, 
instead of treating them as a single fluid. Secondly, we introduced a new scalar force that couples 
differently to CDM and BM. As was discussed recently in 26 , an accurate computation of the 
evolution of the different components is mandatory if one wants to achieve high precision modeling 
of structure formation. In order to test the accuracy of the TRG method for this more general class 
of cosmological scenaria, we analyzed the same cosmologies considered in [40] . and compared our 
results with those of the N-body simulations presented there. The agreement is very good up to 
k ~ 0.5h/Mpc at z = 0, where the non-linear power spectrum is roughly twice the linear one. These 
results confirm the reliability of the TRG as a computational tool that can fill the gap between linear 
and non-linear scales. 

Even in the absence of an extra force, or in the case that the extra force couples universally to all 
matter species, the evolution of BM and CDM differs as a consequence of different initial conditions 
after decoupling. In linear theory, the initial unbalance between BM and CDM fluctuations is almost 
completely washed out by the present epoch. However, when non-linearities are taken into account 
the bias persists. We have seen that the effect is at the percent level in the BAO range at low redshift 
in the uncoupled case, but it may grow up to the 2.5% level when a non-zero coupling is turned on 
with a value compatible with present bounds (obtained within the linear approximation 38, 39 ). 
As a result, these models will receive significant constraints from future galaxy surveys, which aim 
to measure the power spectrum within the BAO range with an accuracy at the percent level. A 
thorough investigation of this model-dependent issue goes beyond the purpose of this paper and is 
postponed for future work. 

The most interesting outcome of our analysis from the point of view of observations concerns 
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Figure 6: The bias in the BAO region for z = 1.12. 



the form of the bias between CDM and BM spectra. In linear theory, the bias is scale indepen- 
dent in recent epochs [25] ■ The non-linear corrections remove this feature, and make the bias scale 
dependent, consistently with the results of [30]. We confirm and extend these results for the BAO 
region (0.03 fo/Mpc Ss k 0.25 fo/Mpc), where we expect the TRG method to be reliable. The effect 
becomes more pronounced with increasing coupling between CDM and DE (assuming that the re- 
spective coupling for BM vanishes). The form of the bias provides an additional observational handle 
for the differentiation between various models. For example, in models with massive neutrinos the 
growth of the spectrum induced by the non-linearities at small length scales is compensated by the 
free-streaming of neutrinos [B]. On the other hand, the bias is expected to remain scale-dependent at 
the non-linear level if there is a substantial CDM-DE coupling. 
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5 Appendix A: Non-linear evolution of perturbations 

We consider first the possibility that the energy density of the Universe is dominated by two com- 
ponents: a) a species of weakly interacting, massive particles, which we identify with dark matter 
(CDM), and b) a slowly varying, classical scalar field <j>, whose contribution is characterized as dark 
energy (DE). There is a direct coupling between the particles and the scalar field: The mass m of the 
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particles depends on the value of (f>. For classical particles, the action of the system can be written as 



with dsi = w g,j, v (xi)dx^ dx\ and the second integral taken over particle trajectories. Variation of the 
action with respect to <j> results in the equation of motion 



dx^ \ d x " J # \Z~5 

with the energy-momentum tensor associated with the gas of particles given by 

{TcdmY" = ^=Y. I ds * ™(<K*i)) ^^8 {A \x~Xi). (5.3) 

Y Q , J CLSi uSi 

Here M is the reduced Planck mass M = (87rG)" 1/2 and we have defined P{cf>)/M = -dlnm(0)/#. 
In the following we shall use eq. (|5.2|l . but we shall approximate the energy- momentum tensor as 
that of an ideal pressureless fluid. The same equation is obtained for scalar-tensor theories of gravity 
in the Einstein frame. If more than one massive species are present in such theories, their coupling 
to the scalar field is universal. In order to be as general as possible when we discuss multiple massive 
species in the following appendix, we assume that each species has a different coupling /3i, induced 
through the ^-dependence of its mass. 

For the metric, we consider an ansatz of the form (|2.2|l . We assume that the Newtonian potential 
$ is weak, $ < 1, and that the field <j> can be decomposed as in eq. (|2.3[) . with 5cj>/(j> <C 1. In 
general, (f> = O(M). The magnitude of the fluctuations of <j) is expected to be similar to that of the 
gravitational field $. The reason is that the source for both is the dark matter density, to which they 
couple with comparable strength. Finally, the density can be decomposed as in eq. (|2.4p . Our aim is 
to take into account the effect of the local velocity field Sv, when this becomes significant because of 
large field gradients. For subhorizon perturbations with momenta k 2> H = a/ a, the linear analysis 
predicts Sv ~ (fc/H)<3? ~ (H/k)(5p/p). We assume that these relations are approximately valid 
even at the non-linear level, within the range of validity of our scheme. Our assumptions can be 
summarized in the hierarchy of scales: $,50/0 <C \Sv\ Sp/p < 1. At the linear level, we have 
Sv 2 ~ &(6p/p). We assume that such a relation holds at the non-linear level as well. As we are 
dealing with subhorizon perturbations, we expect that the spatial derivatives of 5(f) dominate over 
their time derivatives. Following linear theory, we make the more specific assumption that a spatial 
derivative acting on 8(f) or Sv increases the position of that quantity in the hierarchy by one level. 
In this sense V$ is comparable to Sv, while V 2 $ is comparable to p. 

We approximate the energy-momentum tensor of dark matter as (Tcdm)^ v = pVV. We define 
the peculiar velocity through V — Sv" 1 /a. Keeping the leading terms in our expansion gives 

(T C dm)° = p(l + Sv 2 ) 
{Tcdm) 1 \ = —pSv 1 

{TcDuYj = -~p8v l 8vj, (5.4) 

with Svj = Sv 3 . We emphasize at this point that our assumption for the form of (Tcdm)* 3 is consistent 
with the presence of only one gravitational potential $ in our ansatz (|2.2[) for the metric, within the 
order that this potential will be determined through the Einstein equations. 
The energy-momentum tensor of the scalar field has the leading part 



2a 2 

(Ts)\ = 
(Tsf, 



2a 2 



F + U{4>) 



5%, (5.5) 
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with a dot denoting a derivative with respect to r. It also includes a perturbation 



mr = - 2 



. 1^-2 , 2 ^(0) 



i 1^-2 , 2 dim 



(5.6) 



Notice that we have included a term (V50) 2 /2, as it is comparable to Sv 2 or $ within our assumed 
hierarchy. The same is not true for (<50) 2 /2, which is subleading. 

The equation of motion of the scalar field (|5.2|l can be split into an equation for the homogeneous 
part: 

(5-7) 



a am M 



and one for the perturbation: 



i + 2^ 



(/ d0 2 acp M M d<p 



a 2 (p + 5p) 5<f). (5.8) 



The terms with time derivatives are subdominant according to our assumptions. Neglecting them 
results in 



' + 



2 d 2 [/(0) l d/3(0) 



d0 2 



A I 



(P + Sp) 



(5.9) 



For the scalar field to evolve at cosmological times, it must have a mass term d 2 (7(0)/d0 2 = 0(H 2 ). 
We also have (7(0), p = 0(H 2 M 2 ). It is natural to expect df/(0)/d0 = 0(H 2 M). For subhorizon 
perturbations with momenta k 3> H, neglecting the subleading terms results in a very simple Poisson 
equation for the field <50: 



V 2 S<h = 



M H 



(5.10) 



The equation of motion for the gravitational potential $ can be obtained from the first Einstein 
equation. The leading terms give 



3M 2 



a p + -<j) + a (7(0) 



(5.11) 



while the equation for the perturbation is 
V 2 $ - - 3H 2 $ ' 



2M 2 



a 2 Sp + a 2 p 5v 2 — 2< I> + 



+ a 



dU(<t>) 



(5.12) 



Our assumptions about the hierarchy of the various scales and the dominance of the spatial derivatives 
lead to the Poisson equation for the gravitational field $: 



V 2 $ 



2M 2 



— a Sp. 



(5.13) 



We now turn to equations derived from the conservation of the total energy-momentum tensor 
T^ v v = 0. For p = 0, the leading terms give 



p + 3Hp = —^ 0+2H0+a 



d</> J 



M 



(5.14) 



where we have employed eq. (|5.7p . The equation for the perturbations is more complicated. It can be 
simplified considerably through our assumptions about the hierarchy of the various fields. We obtain 



Sp + mSp + V [{p + 8p)5v\ = V 2 <50 = 



M v 



5p, 



(5.15) 
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where we have employed eq. ()5. 10p . 

For /i = i = 1,2, 3 we obtain the generalization of the Euler equation for this system. After 
eliminating higher-order terms we find 

{p + Sp ) [si)+"Hdv+ (Svv) Sv] + ^p + 3Hp + 5p + 3H8p + V[{p + 5p)5v^5v 

= -(p + <5p)V$ + ^ ($ + 2U$- v2 ^ + fl2 ^^) V54>. (5.16) 
Employing eqs. 1(577)1, ((57117)) . (JSTJJ}, (PH5)l we find 

<W + (w-^^0J<W + (<WV) <W=-V$ + ^V<5(/). (5.17) 



For constant /?, we can form a combination of eqs. (|5.10|) . (|5.13)l of the form 



V 2 $ = -^a 2 8p, (5.18) 
2M 

where <3> = $ - (35(f>/M and Af 2 = M 2 /(l + 2/3 2 ). Eq. ()57T7)) becomes 

+ (n- Sv+ (Sv- v) Sv = -V$. (5.19) 

We can see that the Newtonian potential for the perturbations involves a stronger Newton's constant 
G — (87rAf 2 ) -1 = (1 + 2/3 2 )G. There is also a correction ~ cj>Sv in the l.h.s. of the Euler equation, 
because the particles do not follow geodesic motion with respect to the background metric. 



6 Appendix B: Several species of non-relativistic matter 



In this appendix we generalize the formalism to the case that there are more than one particle species 
contributing significantly to the energy density. Each of them couples to the scalar field with a 
different coupling /3i. One example of particular interest includes CDM with /3c dm and baryonic 
matter (BM), for which we assume that Pbm — in order to be consistent with observational 
constraints. We also normalize all dimensionful quantities in terms of the Planck mass. This is 
equivalent to setting M = 1 in the expressions of appendix A. 
The background equations are 



4> + 2H4> + a 



sdU{4>) 



J2p. + U 2 + a 2 U(& 



Pi + 3Hpi = —(f) 'Pi Pi 



(6.1) 

(6.2) 
(6.3) 



with the index i counting the various species. 
For the perturbations, eq. H5.1U|) becomes 



'^2/3iSp 1 , 



(6.4) 



and eq. (|5.13|) 

Eq. (|5.15p generalizes to 



V 2 <£> = ±a 2 ^5 Pl 



(6.5) 



5p\ + 3H5pi + v [{pi + 5pi)8vi] = —0i4>6pi. 



(6.6) 
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Finally, eq. (|5.17p becomes 

Siti + { U - fc(f> ) Svi + { Svi ■ V ) Sv t = - V$ + /3iV5(/>. 



(6.7) 



We can write the above equations in a more useful form by defining the density contrasts Si 
5 pi/ pi < f and 8i(k, r) = V ■ Svi(k, r). For the density contrasts we obtain 



5i + V[(l + 5i)6vi\ =0. 



(6.8) 



For the Fourier transformed quantities, eq. (|6.8[) gives eq. (|2.8|l . while eq. (|2.7b|l gives eq. (|2.10|l . 

We replace time by the variable n — In a(r). For n species of non-relativistic matter we define the 
field ip(k, 77) as a vector with 2n components: 



V2(k,7?) 

<^2n-l(k, ?7) 
V V?2n(k,77) / 



fli(M) 
u 



\ 



n 



This allows us to bring eqs. I|2.8|l . (|2.f 0|l in the usual form [TJ [4J [5] : 

d v (p a (k, rf) + fl ab ip b (k, rj) = e I, 7 a6c (k, -ki, -k2)¥> a (ki, rj)tpbQt2, n). 



(6.9) 



The indices a, b, c take values 1, . . . , In. Repeated momenta are integrated over while repeated indices 
are summed over. The functions 7, that determine effective vertices, are analogous to those employed 
in [?] [5] . We find that the non-zero functions are 

72i-i,2i,2i-i(k,ki,k 2 ) = ^%^^fc(k + ki +k 2 ) = 72i-i,2i-i,2i(k,k 2 ,ki) 

2 (6.f0) 

72 l ,2i, 2l (k,ki,k 2 ) = /3(ki,k 2 ) fo(k + ki + k a)- 

The f2-matrix of eq. (|6.9[) is a 2n x 2n matrix. Let us define the 2x2 matrices 0^(77) and uji^{rj), 
with i, j — 1, . . . , n and i 7^ j, as 



0Ji{v) = 



-1 



and 











(6.11) 



(6.12) 



y --^(2/3^ + 1) 
where a prime denotes a derivative with respect to 77. Then, the fi-matrix can be written as 

( OJ! Wi,2 . . . Wl,n \ 



LJ-> 



L02,: 



V 



(6.13) 



Notice that the only way in which different species of matter influence each other is through the 
matrices oJij, while the vertices do not mix contributions from different species. 
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